% reads in Ds for fractional crystallizaiton

% SBK NEEDS TO FIX THIS SCRIPT. EMAIL ME IF ITS URGENT!
close all 
clear all
addpath('Z_NeededScripts')
addpath('Z_NeededDatasets')
%ColorsHERE
% designate worksheet name to save results:
worksheetName2Save='Ds_FracCryst'; 
worksheetName2Save45='Ds_FracCryst45';
worksheetName2Save='Ds_FracCryst45_HiNiD'; 

%Co_Name = 'PUM'; % This defaults to Hart and Zindler PUM
%Co_Name = 'One'; % This defaults to Hart and Zindler PUM

%%
% Import matrix with mantle bulk compositions
targetStrings_Majors_PETROGEN = {'SiO2','TiO2','Al2O3','Cr2O3','FeO','MgO','CaO','K2O','Na2O'};
DesiredDataName = 'USE'; 
[BulkCompositions, BulkCompName] = reorderElements_Rows('Ds_RevPet', 'MantleBulkCompositions', DesiredDataName, targetStrings_Majors_PETROGEN,'FirstMajor','LastMajor'); 
BulkCompositionsNormalized = sumMgNumNormNoPT(BulkCompositions);
% Import matrix with trace element order and labels
[xlsNumbers, xlsText,xlsRAW] = xlsread('AAA_Step0_PETROGEN_setup_NiD', 'TraceElementOrder');
[ElementRow,ElementColumn]= find(strcmp(xlsRAW,'Elements')==1); 
[LastElementRow,LastElementColumn]= find(strcmp(xlsRAW,'LastElement')==1); 
targetStrings_Trace = xlsRAW(ElementRow+1:LastElementRow-1, ElementColumn)';

U_index = find(strcmp(targetStrings_Trace,'U'));
Th_index = find(strcmp(targetStrings_Trace,'Th'));
Ra_index = find(strcmp(targetStrings_Trace,'Ra'));
H2O_index = find(strcmp(targetStrings_Trace,'H2O'));

% Import mantle source composition. Note that this is not particuarlly important; melt compositions can be  
% converted afterwards to have had a difference source
% [CoRow,CoColumn]= find(strcmp(xlsRAW,Co_Name)==1); 
% Co = xlsRAW(ElementRow+1:LastElementRow-1, CoColumn);
% Co = cell2mat(Co); 

% Import matrix with constant partition coefficients. Note we do change U
% and Th based on pressure
[DRow,DColumn]= find(strcmp(xlsRAW,'FirstD')==1); 
[LastDRow,LastDColumn]= find(strcmp(xlsRAW,'LastD')==1); 
Dmatrix = xlsRAW(ElementRow+1:LastElementRow-1, DColumn:LastDColumn);
Dmatrix = cell2mat(Dmatrix); 

D_Th = eval(sprintf('Dmatrix(%d,:)',Th_index));
D_U = eval(sprintf('Dmatrix(%d,:)',U_index));
D_U./D_Th;
worksheetName2Save

save(worksheetName2Save45,'Dmatrix')

%Dmatrix(45,:)=[]; 
save(worksheetName2Save,'Dmatrix')